library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)library(readr)library(tidyr)# Read the CAPTAIN2 EDGE2 RDS fileCAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates100.rds"))# Analyze non-zero cells in the prioritization outputcat("Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:\n")
Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:
cat("Cells with zero priority:", zero_cells, " (", round(zero_cells/total_cells*100, 2), "%)\n", sep="")
Cells with zero priority:226788 (97.52%)
Code
cat("Cells with NA priority:", na_cells, " (", round(na_cells/total_cells*100, 2), "%)\n", sep="")
Cells with NA priority:0 (0%)
Code
# Summary statistics of priority valuespriority_summary <-summary(CAPTAIN2_EDGE2_data$Priority)cat("\nSummary statistics of priority values:\n")
Summary statistics of priority values:
Code
print(priority_summary)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00000 0.00000 0.00000 0.02138 0.00000 1.00000
Code
# Distribution of non-zero priority valuesnonzero_priority <- CAPTAIN2_EDGE2_data$Priority[CAPTAIN2_EDGE2_data$Priority >0]cat("\nDistribution of non-zero priority values:\n")
# Load one of your input raster files to extract the correct grid structureraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 data based on PUIDCAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_EDGE2_sf <-st_as_sf( CAPTAIN2_EDGE2_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_EDGE2 <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_EDGE2 <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_EDGE2 <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_EDGE2_plot <-ggplot() +geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_EDGE2# Display the plotprint(CAPTAIN2_EDGE2_plot)
Code
# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_EDGE2_priorities_01.png"),plot = CAPTAIN2_EDGE2_plot,width =10,height =6,dpi =300,bg ="white")
FUSE continental 0.1 budget
Code
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 FUSE RDS fileCAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 FUSE data based on PUIDCAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_FUSE_sf <-st_as_sf( CAPTAIN2_FUSE_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_FUSE <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_FUSE <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_FUSE <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_FUSE_plot <-ggplot() +geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_FUSE, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_FUSE, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_FUSE# Display the plotprint(CAPTAIN2_FUSE_plot)
Code
# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_FUSE_priorities_01.png"),plot = CAPTAIN2_FUSE_plot,width =10,height =6,dpi =300,bg ="white")
IUCN continental 0.1 budget
Code
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 IUCN RDS fileCAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 IUCN data based on PUIDCAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_IUCN_sf <-st_as_sf( CAPTAIN2_IUCN_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_IUCN <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_IUCN <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_IUCN <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_IUCN_plot <-ggplot() +geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_IUCN, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_IUCN, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_IUCN# Display the plotprint(CAPTAIN2_IUCN_plot)
Code
# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_priorities_01.png"),plot = CAPTAIN2_IUCN_plot,width =10,height =6,dpi =300,bg ="white")
Difference maps
Code
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read all three index RDS filesCAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates100.rds"))CAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates100.rds"))CAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structureraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Join all three datasets with coordinatesIUCN_with_coords <- CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority) %>%left_join(coords, by ="PUID")EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority) %>%left_join(coords, by ="PUID") FUSE_with_coords <- CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority) %>%left_join(coords, by ="PUID")# Combine all datasetsall_indices <- coords %>%left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by ="PUID") %>%left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by ="PUID") %>%left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by ="PUID")# Calculate differencesall_indices <- all_indices %>%mutate(# Replace NA with 0 for calculation purposesIUCN =ifelse(is.na(IUCN), 0, IUCN),EDGE2 =ifelse(is.na(EDGE2), 0, EDGE2),FUSE =ifelse(is.na(FUSE), 0, FUSE),# Calculate differencesIUCN_minus_FUSE = IUCN - FUSE,IUCN_minus_EDGE2 = IUCN - EDGE2,EDGE2_minus_FUSE = EDGE2 - FUSE )# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected <-st_transform(world, crs = mcbryde_thomas_2)# Create base thememy_theme <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Filter to non-zero differences for each comparison to reduce plot size# IUCN - FUSEIUCN_FUSE_diff <- all_indices %>%filter(IUCN_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# IUCN - EDGE2IUCN_EDGE2_diff <- all_indices %>%filter(IUCN_minus_EDGE2 !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# EDGE2 - FUSEEDGE2_FUSE_diff <- all_indices %>%filter(EDGE2_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# Create a diverging color palette for difference maps# Blue for negative (first index lower), white for zero, red for positive (first index higher)diff_colors <-c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")# 1. IUCN - FUSE Difference MapIUCN_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus FUSE Index",x =NULL, y =NULL) + my_theme# 2. IUCN - EDGE2 Difference MapIUCN_EDGE2_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - EDGE2)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus EDGE2 Index",x =NULL, y =NULL) + my_theme# 3. EDGE2 - FUSE Difference MapEDGE2_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(EDGE2 - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="EDGE2 Index minus FUSE Index",x =NULL, y =NULL) + my_theme# Display all plotsprint(IUCN_FUSE_plot)
Code
print(IUCN_EDGE2_plot)
Code
print(EDGE2_FUSE_plot)
Code
# Save all plotsggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_minus_FUSE_difference.png"),plot = IUCN_FUSE_plot,width =10,height =6,dpi =300,bg ="white")ggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_minus_EDGE2_difference.png"),plot = IUCN_EDGE2_plot,width =10,height =6,dpi =300,bg ="white")ggsave(filename = here::here("outputs", "CAPTAIN2_EDGE2_minus_FUSE_difference.png"),plot = EDGE2_FUSE_plot,width =10,height =6,dpi =300,bg ="white")# Optionally, create a panel with all three difference mapslibrary(patchwork)# Combine all plotsall_diffs_plot <- IUCN_FUSE_plot / IUCN_EDGE2_plot / EDGE2_FUSE_plot +plot_annotation(title ="Differences Between Conservation Priority Indices",subtitle ="Budget: 0.1, Replicates: 100",theme =theme(plot.title =element_text(hjust =0.5),plot.subtitle =element_text(hjust =0.5)) )#all_diffs_plot# Save the combined plot#ggsave(# filename = here::here("outputs", "CAPTAIN2_all_differences.png"),# plot = all_diffs_plot,# width = 10,# height = 15,# dpi = 300,# bg = "white"#)
Species level priorities
Code
# Load required packageslibrary(here)library(readr)library(dplyr)library(tidyr)library(ggplot2)# Read the protected range fractions RDS fileprotected_fractions <-readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions.rds"))# Read the continental shark conservation metrics CSV fileshark_metrics <-read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))# Rename column in shark_metrics to match bettershark_metrics <- shark_metrics %>%rename(Species =`Species name`)# Order shark_metrics alphabetically by Species nameshark_metrics <- shark_metrics %>%arrange(Species)# Add an original order ID to protected_fractions to maintain its original orderprotected_fractions$original_order <-1:nrow(protected_fractions)# Add row number as IDs to both datasetsprotected_fractions$protected_ID <-1:nrow(protected_fractions)shark_metrics$species_ID <-1:nrow(shark_metrics)# Check if the datasets have the same number of rowsif(nrow(protected_fractions) ==nrow(shark_metrics)) {# Create index mapping - this maintains the original protection data ordering# while allowing us to associate with alphabetically ordered species names indices <-data.frame(protected_ID =1:nrow(protected_fractions),species_ID =1:nrow(shark_metrics) )# Join protected_fractions with indices protected_with_indices <- protected_fractions %>%left_join(indices, by ="protected_ID")# Join shark_metrics with indices shark_with_indices <- shark_metrics %>%left_join(indices, by ="species_ID")# Now join the datasets, matching on species_ID and protected_ID combined_data <- protected_with_indices %>%inner_join( shark_with_indices,by =c("species_ID", "protected_ID"),suffix =c("_captain", "_original") ) %>%# Sort by the original order of protected_fractions arrange(original_order)cat("Successfully joined datasets with", nrow(combined_data), "species\n")cat("First few species in combined dataset:\n")print(head(combined_data[, c("Species_captain", "Species_original")]))# Define IUCN categories and order - using only the first 5 categories iucn_labels <-c("1"="LC", "2"="NT", "3"="VU", "4"="EN", "5"="CR" ) iucn_order <-c("LC", "NT", "VU", "EN", "CR")# Define colors for IUCN categories iucn_colors <-c("LC"="#50C878", # Green"NT"="#FFFF00", # Yellow"VU"="#FFA500", # Orange"EN"="#FF8C00", # Dark Orange"CR"="#FF0000"# Red )# 1. IUCN Boxplot iucn_boxplot <- combined_data %>%mutate(IUCN_status =factor(IUCN_original, levels =1:5, labels = iucn_order),protection_percentage = IUCN_captain *100 ) %>%ggplot(aes(x = IUCN_status, y = protection_percentage)) +# geom_violin(aes(fill = IUCN_status, color = IUCN_status), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="IUCN Priority Index: Range Protection by IUCN Status",x ="IUCN Red List threat status", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_fill_manual(values = iucn_colors) +scale_color_manual(values = iucn_colors) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 2. FUSE Boxplot fuse_boxplot <- combined_data %>%mutate(FUSE_category =factor(FUSE_original),protection_percentage = FUSE_captain *100 ) %>%ggplot(aes(x = FUSE_category, y = protection_percentage)) +# geom_violin(aes(fill = FUSE_category, color = FUSE_category), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="FUSE Priority Index: Range Protection by FUSE Score",x ="FUSE Score", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 3. EDGE2 Boxplot edge2_boxplot <- combined_data %>%mutate(EDGE2_category =factor(EDGE2_original),protection_percentage = EDGE2_captain *100 ) %>%ggplot(aes(x = EDGE2_category, y = protection_percentage)) +# geom_violin(aes(fill = EDGE2_category, color = EDGE2_category), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="EDGE2 Priority Index: Range Protection by EDGE2 Score",x ="EDGE2 Score", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# Display individual plotsprint(iucn_boxplot)print(fuse_boxplot)print(edge2_boxplot)# Save individual plotsggsave(here::here("outputs", "iucn_protection_boxplot.png"), iucn_boxplot, width =10, height =6, dpi =300)ggsave(here::here("outputs", "fuse_protection_boxplot.png"), fuse_boxplot, width =10, height =6, dpi =300)ggsave(here::here("outputs", "edge2_protection_boxplot.png"), edge2_boxplot, width =10, height =6, dpi =300)# Create a nicely formatted table showing the ordered species species_table <- combined_data %>% dplyr::select(Species_captain, Species_original, IUCN_captain, FUSE_captain, EDGE2_captain, IUCN_original, FUSE_original, EDGE2_original) %>%arrange(Species_original) %>%head(20) # Just show the first 20 for display# Print species tablecat("\nFirst 20 species (alphabetically by original species name):\n")print(species_table)# Save the full combined datawrite.csv(combined_data, here::here("outputs", "combined_species_data.csv"), row.names =FALSE)} else {cat("ERROR: Datasets have different number of rows.\n")cat("Protected fractions:", nrow(protected_fractions), "rows\n")cat("Shark metrics:", nrow(shark_metrics), "rows\n")}
Successfully joined datasets with 1000 species
First few species in combined dataset:
Species_captain Species_original
1 Species_1 Acroteriobatus_annulatus
2 Species_2 Acroteriobatus_blochii
3 Species_3 Acroteriobatus_leucospilus
4 Species_4 Acroteriobatus_ocellatus
5 Species_5 Acroteriobatus_salalah
6 Species_6 Acroteriobatus_variegatus
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 FUSE RDS fileCAPTAIN2_res_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 FUSE data based on PUIDCAPTAIN2_res_FUSE_data_with_coords <- CAPTAIN2_res_FUSE_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_res_FUSE_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_res_FUSE_data_nonzero <- CAPTAIN2_res_FUSE_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_res_FUSE_sf <-st_as_sf( CAPTAIN2_res_FUSE_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_res_FUSE <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_res_FUSE <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_res_FUSE <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_res_FUSE_plot <-ggplot() +geom_sf(data = CAPTAIN2_res_FUSE_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_res_FUSE, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_res_FUSE, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_res_FUSE# Display the plotprint(CAPTAIN2_res_FUSE_plot)
Code
# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_res_FUSE_priorities_01.png"),plot = CAPTAIN2_res_FUSE_plot,width =10,height =6,dpi =300,bg ="white")
Species level priorities
Code
# Load required packageslibrary(here)library(readr)library(dplyr)library(ggplot2)# Read the protected range fractions RDS file (with "FUSE_res")protected_fractions <-readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_with_FUSE_res.rds"))# Read the continental shark conservation metrics CSV fileshark_metrics <-read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))# Rename column in shark_metrics to match bettershark_metrics <- shark_metrics %>%rename(Species =`Species name`)# Order shark_metrics alphabetically by Species nameshark_metrics <- shark_metrics %>%arrange(Species)# Add an original order ID to protected_fractions to maintain its original orderprotected_fractions$original_order <-1:nrow(protected_fractions)# Add row number as IDs to both datasetsprotected_fractions$protected_ID <-1:nrow(protected_fractions)shark_metrics$species_ID <-1:nrow(shark_metrics)# Check if the datasets have the same number of rowsif (nrow(protected_fractions) ==nrow(shark_metrics)) {# Create index mapping - this maintains the original protection data ordering# while allowing us to associate with alphabetically ordered species names indices <-data.frame(protected_ID =1:nrow(protected_fractions),species_ID =1:nrow(shark_metrics) )# Join protected_fractions with indices protected_with_indices <- protected_fractions %>%left_join(indices, by ="protected_ID")# Join shark_metrics with indices shark_with_indices <- shark_metrics %>%left_join(indices, by ="species_ID")# Now join the datasets, matching on species_ID and protected_ID combined_data <- protected_with_indices %>%inner_join( shark_with_indices,by =c("species_ID", "protected_ID"),suffix =c("_captain", "_original") ) %>%# Sort by the original order of protected_fractions arrange(original_order)cat("Successfully joined datasets with", nrow(combined_data), "species\n")# Create the FUSE_res vs FUSE_original Boxplot fuse_res_vs_fuse_original_plot <- combined_data %>%mutate(FUSE_category =factor(FUSE_original), # Convert FUSE_original to a factor for x-axis categoriesprotection_percentage = FUSE_res *100# Convert FUSE_res to percentage ) %>%ggplot(aes(x = FUSE_category, y = protection_percentage)) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.3, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="FUSE_res Priority Index: Range Protection by FUSE Category",x ="FUSE Original Category", y ="Range Protected (%)" ) +theme_minimal() +theme(plot.title =element_text(hjust =0.5),panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25) )# Display the plotprint(fuse_res_vs_fuse_original_plot)# Save the plotggsave(filename = here::here("outputs", "fuse_res_vs_fuse_original_boxplot.png"),plot = fuse_res_vs_fuse_original_plot,width =10,height =6,dpi =300 )# Save the full combined datawrite.csv( combined_data, here::here("outputs", "combined_species_data_with_FUSE_res.csv"), row.names =FALSE )cat("\nFUSE_res vs FUSE_original plot saved as 'fuse_res_vs_fuse_original_boxplot.png'\n")cat("Combined dataset saved as 'combined_species_data_with_FUSE_res.csv'\n")} else {cat("ERROR: Datasets have different number of rows.\n")cat("Protected fractions:", nrow(protected_fractions), "rows\n")cat("Shark metrics:", nrow(shark_metrics), "rows\n")}
Successfully joined datasets with 1000 species
FUSE_res vs FUSE_original plot saved as 'fuse_res_vs_fuse_original_boxplot.png'
Combined dataset saved as 'combined_species_data_with_FUSE_res.csv'
Source Code
---title: "CAPTAIN2 workflow"author: "Théophile L. Mouton"date: "May 21, 2025"format: html: toc: true toc-location: right css: custom.css output-file: "CAPTAIN2 workflow.html" self-contained: true code-fold: true code-tools: trueeditor: visualexecute: warning: false message: false echo: true---## EDGE2 continental 0.1 budget```{r}library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)library(readr)library(tidyr)# Read the CAPTAIN2 EDGE2 RDS fileCAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates100.rds"))# Analyze non-zero cells in the prioritization outputcat("Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:\n")total_cells <-nrow(CAPTAIN2_EDGE2_data)nonzero_cells <-sum(CAPTAIN2_EDGE2_data$Priority >0, na.rm =TRUE)zero_cells <-sum(CAPTAIN2_EDGE2_data$Priority ==0, na.rm =TRUE)na_cells <-sum(is.na(CAPTAIN2_EDGE2_data$Priority))cat("Total cells in grid:", total_cells, "\n")cat("Cells with non-zero priority:", nonzero_cells, " (", round(nonzero_cells/total_cells*100, 2), "%)\n", sep="")cat("Cells with zero priority:", zero_cells, " (", round(zero_cells/total_cells*100, 2), "%)\n", sep="")cat("Cells with NA priority:", na_cells, " (", round(na_cells/total_cells*100, 2), "%)\n", sep="")# Summary statistics of priority valuespriority_summary <-summary(CAPTAIN2_EDGE2_data$Priority)cat("\nSummary statistics of priority values:\n")print(priority_summary)# Distribution of non-zero priority valuesnonzero_priority <- CAPTAIN2_EDGE2_data$Priority[CAPTAIN2_EDGE2_data$Priority >0]cat("\nDistribution of non-zero priority values:\n")priority_quantiles <-quantile(nonzero_priority, probs =seq(0, 1, 0.1), na.rm =TRUE)print(priority_quantiles)# Load one of your input raster files to extract the correct grid structureraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 data based on PUIDCAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_EDGE2_sf <-st_as_sf( CAPTAIN2_EDGE2_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_EDGE2 <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_EDGE2 <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_EDGE2 <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_EDGE2_plot <-ggplot() +geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_EDGE2# Display the plotprint(CAPTAIN2_EDGE2_plot)# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_EDGE2_priorities_01.png"),plot = CAPTAIN2_EDGE2_plot,width =10,height =6,dpi =300,bg ="white")```# FUSE continental 0.1 budget```{r}library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 FUSE RDS fileCAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 FUSE data based on PUIDCAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_FUSE_sf <-st_as_sf( CAPTAIN2_FUSE_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_FUSE <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_FUSE <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_FUSE <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_FUSE_plot <-ggplot() +geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_FUSE, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_FUSE, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_FUSE# Display the plotprint(CAPTAIN2_FUSE_plot)# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_FUSE_priorities_01.png"),plot = CAPTAIN2_FUSE_plot,width =10,height =6,dpi =300,bg ="white")```# IUCN continental 0.1 budget ```{r}library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 IUCN RDS fileCAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 IUCN data based on PUIDCAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_IUCN_sf <-st_as_sf( CAPTAIN2_IUCN_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_IUCN <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_IUCN <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_IUCN <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_IUCN_plot <-ggplot() +geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_IUCN, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_IUCN, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_IUCN# Display the plotprint(CAPTAIN2_IUCN_plot)# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_priorities_01.png"),plot = CAPTAIN2_IUCN_plot,width =10,height =6,dpi =300,bg ="white")```# Difference maps ```{r}library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read all three index RDS filesCAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates100.rds"))CAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates100.rds"))CAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structureraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Join all three datasets with coordinatesIUCN_with_coords <- CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority) %>%left_join(coords, by ="PUID")EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority) %>%left_join(coords, by ="PUID") FUSE_with_coords <- CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority) %>%left_join(coords, by ="PUID")# Combine all datasetsall_indices <- coords %>%left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by ="PUID") %>%left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by ="PUID") %>%left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by ="PUID")# Calculate differencesall_indices <- all_indices %>%mutate(# Replace NA with 0 for calculation purposesIUCN =ifelse(is.na(IUCN), 0, IUCN),EDGE2 =ifelse(is.na(EDGE2), 0, EDGE2),FUSE =ifelse(is.na(FUSE), 0, FUSE),# Calculate differencesIUCN_minus_FUSE = IUCN - FUSE,IUCN_minus_EDGE2 = IUCN - EDGE2,EDGE2_minus_FUSE = EDGE2 - FUSE )# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected <-st_transform(world, crs = mcbryde_thomas_2)# Create base thememy_theme <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Filter to non-zero differences for each comparison to reduce plot size# IUCN - FUSEIUCN_FUSE_diff <- all_indices %>%filter(IUCN_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# IUCN - EDGE2IUCN_EDGE2_diff <- all_indices %>%filter(IUCN_minus_EDGE2 !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# EDGE2 - FUSEEDGE2_FUSE_diff <- all_indices %>%filter(EDGE2_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# Create a diverging color palette for difference maps# Blue for negative (first index lower), white for zero, red for positive (first index higher)diff_colors <-c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")# 1. IUCN - FUSE Difference MapIUCN_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus FUSE Index",x =NULL, y =NULL) + my_theme# 2. IUCN - EDGE2 Difference MapIUCN_EDGE2_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - EDGE2)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus EDGE2 Index",x =NULL, y =NULL) + my_theme# 3. EDGE2 - FUSE Difference MapEDGE2_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(EDGE2 - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="EDGE2 Index minus FUSE Index",x =NULL, y =NULL) + my_theme# Display all plotsprint(IUCN_FUSE_plot)print(IUCN_EDGE2_plot)print(EDGE2_FUSE_plot)# Save all plotsggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_minus_FUSE_difference.png"),plot = IUCN_FUSE_plot,width =10,height =6,dpi =300,bg ="white")ggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_minus_EDGE2_difference.png"),plot = IUCN_EDGE2_plot,width =10,height =6,dpi =300,bg ="white")ggsave(filename = here::here("outputs", "CAPTAIN2_EDGE2_minus_FUSE_difference.png"),plot = EDGE2_FUSE_plot,width =10,height =6,dpi =300,bg ="white")# Optionally, create a panel with all three difference mapslibrary(patchwork)# Combine all plotsall_diffs_plot <- IUCN_FUSE_plot / IUCN_EDGE2_plot / EDGE2_FUSE_plot +plot_annotation(title ="Differences Between Conservation Priority Indices",subtitle ="Budget: 0.1, Replicates: 100",theme =theme(plot.title =element_text(hjust =0.5),plot.subtitle =element_text(hjust =0.5)) )#all_diffs_plot# Save the combined plot#ggsave(# filename = here::here("outputs", "CAPTAIN2_all_differences.png"),# plot = all_diffs_plot,# width = 10,# height = 15,# dpi = 300,# bg = "white"#)```# Species level priorities ```{r}# Load required packageslibrary(here)library(readr)library(dplyr)library(tidyr)library(ggplot2)# Read the protected range fractions RDS fileprotected_fractions <-readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions.rds"))# Read the continental shark conservation metrics CSV fileshark_metrics <-read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))# Rename column in shark_metrics to match bettershark_metrics <- shark_metrics %>%rename(Species =`Species name`)# Order shark_metrics alphabetically by Species nameshark_metrics <- shark_metrics %>%arrange(Species)# Add an original order ID to protected_fractions to maintain its original orderprotected_fractions$original_order <-1:nrow(protected_fractions)# Add row number as IDs to both datasetsprotected_fractions$protected_ID <-1:nrow(protected_fractions)shark_metrics$species_ID <-1:nrow(shark_metrics)# Check if the datasets have the same number of rowsif(nrow(protected_fractions) ==nrow(shark_metrics)) {# Create index mapping - this maintains the original protection data ordering# while allowing us to associate with alphabetically ordered species names indices <-data.frame(protected_ID =1:nrow(protected_fractions),species_ID =1:nrow(shark_metrics) )# Join protected_fractions with indices protected_with_indices <- protected_fractions %>%left_join(indices, by ="protected_ID")# Join shark_metrics with indices shark_with_indices <- shark_metrics %>%left_join(indices, by ="species_ID")# Now join the datasets, matching on species_ID and protected_ID combined_data <- protected_with_indices %>%inner_join( shark_with_indices,by =c("species_ID", "protected_ID"),suffix =c("_captain", "_original") ) %>%# Sort by the original order of protected_fractions arrange(original_order)cat("Successfully joined datasets with", nrow(combined_data), "species\n")cat("First few species in combined dataset:\n")print(head(combined_data[, c("Species_captain", "Species_original")]))# Define IUCN categories and order - using only the first 5 categories iucn_labels <-c("1"="LC", "2"="NT", "3"="VU", "4"="EN", "5"="CR" ) iucn_order <-c("LC", "NT", "VU", "EN", "CR")# Define colors for IUCN categories iucn_colors <-c("LC"="#50C878", # Green"NT"="#FFFF00", # Yellow"VU"="#FFA500", # Orange"EN"="#FF8C00", # Dark Orange"CR"="#FF0000"# Red )# 1. IUCN Boxplot iucn_boxplot <- combined_data %>%mutate(IUCN_status =factor(IUCN_original, levels =1:5, labels = iucn_order),protection_percentage = IUCN_captain *100 ) %>%ggplot(aes(x = IUCN_status, y = protection_percentage)) +# geom_violin(aes(fill = IUCN_status, color = IUCN_status), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="IUCN Priority Index: Range Protection by IUCN Status",x ="IUCN Red List threat status", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_fill_manual(values = iucn_colors) +scale_color_manual(values = iucn_colors) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 2. FUSE Boxplot fuse_boxplot <- combined_data %>%mutate(FUSE_category =factor(FUSE_original),protection_percentage = FUSE_captain *100 ) %>%ggplot(aes(x = FUSE_category, y = protection_percentage)) +# geom_violin(aes(fill = FUSE_category, color = FUSE_category), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="FUSE Priority Index: Range Protection by FUSE Score",x ="FUSE Score", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 3. EDGE2 Boxplot edge2_boxplot <- combined_data %>%mutate(EDGE2_category =factor(EDGE2_original),protection_percentage = EDGE2_captain *100 ) %>%ggplot(aes(x = EDGE2_category, y = protection_percentage)) +# geom_violin(aes(fill = EDGE2_category, color = EDGE2_category), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="EDGE2 Priority Index: Range Protection by EDGE2 Score",x ="EDGE2 Score", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# Display individual plotsprint(iucn_boxplot)print(fuse_boxplot)print(edge2_boxplot)# Save individual plotsggsave(here::here("outputs", "iucn_protection_boxplot.png"), iucn_boxplot, width =10, height =6, dpi =300)ggsave(here::here("outputs", "fuse_protection_boxplot.png"), fuse_boxplot, width =10, height =6, dpi =300)ggsave(here::here("outputs", "edge2_protection_boxplot.png"), edge2_boxplot, width =10, height =6, dpi =300)# Create a nicely formatted table showing the ordered species species_table <- combined_data %>% dplyr::select(Species_captain, Species_original, IUCN_captain, FUSE_captain, EDGE2_captain, IUCN_original, FUSE_original, EDGE2_original) %>%arrange(Species_original) %>%head(20) # Just show the first 20 for display# Print species tablecat("\nFirst 20 species (alphabetically by original species name):\n")print(species_table)# Save the full combined datawrite.csv(combined_data, here::here("outputs", "combined_species_data.csv"), row.names =FALSE)} else {cat("ERROR: Datasets have different number of rows.\n")cat("Protected fractions:", nrow(protected_fractions), "rows\n")cat("Shark metrics:", nrow(shark_metrics), "rows\n")}```# RES models ## FUSE continental 0.1 budget```{r}library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 FUSE RDS fileCAPTAIN2_res_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates100.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 FUSE data based on PUIDCAPTAIN2_res_FUSE_data_with_coords <- CAPTAIN2_res_FUSE_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_res_FUSE_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_res_FUSE_data_nonzero <- CAPTAIN2_res_FUSE_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_res_FUSE_sf <-st_as_sf( CAPTAIN2_res_FUSE_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_res_FUSE <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_res_FUSE <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_res_FUSE <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_res_FUSE_plot <-ggplot() +geom_sf(data = CAPTAIN2_res_FUSE_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_res_FUSE, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_res_FUSE, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 100",x =NULL, y =NULL) + my_theme_CAPTAIN2_res_FUSE# Display the plotprint(CAPTAIN2_res_FUSE_plot)# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_res_FUSE_priorities_01.png"),plot = CAPTAIN2_res_FUSE_plot,width =10,height =6,dpi =300,bg ="white")```## Species level priorities```{r}# Load required packageslibrary(here)library(readr)library(dplyr)library(ggplot2)# Read the protected range fractions RDS file (with "FUSE_res")protected_fractions <-readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_with_FUSE_res.rds"))# Read the continental shark conservation metrics CSV fileshark_metrics <-read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))# Rename column in shark_metrics to match bettershark_metrics <- shark_metrics %>%rename(Species =`Species name`)# Order shark_metrics alphabetically by Species nameshark_metrics <- shark_metrics %>%arrange(Species)# Add an original order ID to protected_fractions to maintain its original orderprotected_fractions$original_order <-1:nrow(protected_fractions)# Add row number as IDs to both datasetsprotected_fractions$protected_ID <-1:nrow(protected_fractions)shark_metrics$species_ID <-1:nrow(shark_metrics)# Check if the datasets have the same number of rowsif (nrow(protected_fractions) ==nrow(shark_metrics)) {# Create index mapping - this maintains the original protection data ordering# while allowing us to associate with alphabetically ordered species names indices <-data.frame(protected_ID =1:nrow(protected_fractions),species_ID =1:nrow(shark_metrics) )# Join protected_fractions with indices protected_with_indices <- protected_fractions %>%left_join(indices, by ="protected_ID")# Join shark_metrics with indices shark_with_indices <- shark_metrics %>%left_join(indices, by ="species_ID")# Now join the datasets, matching on species_ID and protected_ID combined_data <- protected_with_indices %>%inner_join( shark_with_indices,by =c("species_ID", "protected_ID"),suffix =c("_captain", "_original") ) %>%# Sort by the original order of protected_fractions arrange(original_order)cat("Successfully joined datasets with", nrow(combined_data), "species\n")# Create the FUSE_res vs FUSE_original Boxplot fuse_res_vs_fuse_original_plot <- combined_data %>%mutate(FUSE_category =factor(FUSE_original), # Convert FUSE_original to a factor for x-axis categoriesprotection_percentage = FUSE_res *100# Convert FUSE_res to percentage ) %>%ggplot(aes(x = FUSE_category, y = protection_percentage)) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.3, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="FUSE_res Priority Index: Range Protection by FUSE Category",x ="FUSE Original Category", y ="Range Protected (%)" ) +theme_minimal() +theme(plot.title =element_text(hjust =0.5),panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25) )# Display the plotprint(fuse_res_vs_fuse_original_plot)# Save the plotggsave(filename = here::here("outputs", "fuse_res_vs_fuse_original_boxplot.png"),plot = fuse_res_vs_fuse_original_plot,width =10,height =6,dpi =300 )# Save the full combined datawrite.csv( combined_data, here::here("outputs", "combined_species_data_with_FUSE_res.csv"), row.names =FALSE )cat("\nFUSE_res vs FUSE_original plot saved as 'fuse_res_vs_fuse_original_boxplot.png'\n")cat("Combined dataset saved as 'combined_species_data_with_FUSE_res.csv'\n")} else {cat("ERROR: Datasets have different number of rows.\n")cat("Protected fractions:", nrow(protected_fractions), "rows\n")cat("Shark metrics:", nrow(shark_metrics), "rows\n")}```